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During this six month period of the program investigations were carried out in 
two general areas, (i) Energy barriers and structural transitions between isomers 
of small A1 clusters were investigated. In this study an empirical potential function 
which was parametrized based on accurate first principle results was employed, 
(ii) A comparative study was conducted to investigate the applicability of most 
commonly employed model potential functions in calculating various bulk, surface 
and small cluster properties. 


I. Calculation of Energy Barriers and Transition Structures for Small A1 
Clusters. 

A model potential energy function fitted to high level ab initio results is used 
for the first time to calculate saddle point configurations and energies between 
low lying structures of Al 4 , AI 5 and Alg. For these clusters, energy barriers for 
structural isomerizations were found to be of varying heights. In general, barriers 
between high-symmetry equilibrium structures are low indicating that for those 
cases isomerization reactions in the gas-phase may take place at lower temperatures. 
On the other hand, calculations indicate that barriers between two- and three- 
dimensional forms are high which may affect the equilibration process in the gas- 
phase. 

This work has been presented in the Fourth International Symposium on Small 
Particles and Inorganic Clusters held in Aix-en-Province, France, on July 5-9, 1988. 
Appendix I includes the manuscript which will be published in “Atoms, Molecules 
and Clusters” (Zeitschrift fur Pyhsik D). 
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II. Comparison of Interatomic Potentials Based on Multi-Body Interac- 
tions 

In this study, four most commonly used model potentials comprising two- 
and three-body interactions were analyzed comparatively. We included “Tersoff”, 
“Tersoff-Dodson”, “Stillenger- Weber” and “Lennard- Jones- Axilrod- Teller” poten- 
tials, and investigated their capabilities in reproducing various bulk, surface and 
small cluster properties of some selected elements. Applicability and limitations of 
each potential function were delineated and discussed separately. In general, struc- 
tural properties calculated by these functions were found to be in better agreement 
with experiment than corresponding energy-related properties. 

The first part of this study dealing with properties of silicon has been now 
completed and the manuscript has been sent for publication. A shortend version of 
this manuscript is presented in Appendix II. 
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APPENDIX I. 


Energy Barriers and Structural Transitions of Small A1 Clusters 


Recent experimental and theoretical investigations on small clusters indicate 
that certain microclusters exhibit isomeric transformations in the gas phase [1-7]. 
In order to thoroughly understand and analyze those isomerization processes, all 
the stationary points (at least for the low-lying configurations) of the potential en- 
ergy surface associated with the cluster, must be estimated [8]. Stationary points 
include not only the minima of the surface, but also saddle point configurations 
[9]. At the present time, highly accurate calculations based on ab initio methods, 
are possible only for geometry optimizations of relatively smaller clusters [10]. A 
high level calculation for the lowest saddle point of paths connecting two minima, 
however, is a formidable job even with today’s supercomputers. In this work, using 
a model potential energy function, which takes into account two- and three-body 
interactions, we calculated saddle point structures and corresponding energies be- 
tween several low-lying equilibrium configurations of AI4, Al 5 and Ale- The model 
potential was parametrized by fitting it to highly accurate ab initio results and it 
has been shown that it is able to reproduce many equilibrium structures of small 
A1 clusters (up to 6 atoms) with an acceptable accuracy [10]. The functional form 
of the potential and values of parameters are given in reference 10. Configura- 
tions corresponding to energy minima were calculated in this work, by an energy 
gradient method [9]. Low-lying configurations were selected from a large number 
of results calculated with different initial geometries. Saddle points like minima 
were also treated as stationary points on the energy surface. Several different paths 
between two configurational minima were taken into consideration and the saddle 
point structure corresponding to the path with the lowest barrier is reported in this 
study. Calculated results are summarized in Table 1. 

For the AI4 case, the two lowest minimum energy structures corresponding 
to three-dimensional (3d) and two-dimensional (2d) configurations, are the regular 
tetrahedron (T)4 and the rhombus (R ) 4 structures, respectively. The saddle point 
structure for the (T ) 4 — ♦ (i?)4 transformation is a highly distorted tetrahedron with 
one elongated side. The barrier for this 3d— >2d transition was estimated as 0.0299 
eV /atom. 

Energetically low-lying 3d configurations for the AI5 case are the triangular 
bipyramid ( TBP) 5 which is the most compact structure and the square pyramid 
(SP) 5. The energetically most favorable Al 5 structure is the 2d close-packed (CP ) 5 
configuration. Two saddle point structures were calculated, in this case. For the 
transition between two 3d structures, the saddle point configuration was found to 
be a distorted triangular bipyramid. Its two apex atoms are shifted toward to one 
side of the triangular base. Calculated energy for this saddle point configuration 
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corresponds to a 0.0055 eV/atom barrier height for the (TBP) 5 — ► (SP) 5 trans- 
formation. For the 3d— >2d transition, among many possible alternatives, the path 
with the lowest saddle point energy was found for the (SP) 5 — ► (CP) 5 transfor- 
mation. This saddle point structure resembles very much the ( CP)$ but it is not 
planar. The barrier for this (SP ) 5 — ► ( CP)s transition was calculated to be 0.0259 
eV/atom. 

For Alg, we considered three low-lying 3d configurations. The square bipyra- 
mid (SBP) 6 structure is the most compact configuration. The tripyramid ( TP) 6 is 
formed by three slightly distorted tetrahedra. Structurally it is similar to (TBP) 5 
with an additional A1 atom positioned above one of its faces. The pentagonal pyra- 
mid (PP) 6 is the energetically most stable 3d configuration. The low-lying 2d con- 
figurations for Al 6 include the close-packed (CP)& and the triangular close-packed 
(TCP) 6 structures. Three transition structures were investigated in this case. The 
saddle point configuration for the (SBP)t — ► (TP)$ transition, resembles the (TP)e 
structure with three distorted tetrahedra. The commonly shared edge of three tetra- 
hedra, in this configuration, was found to be considerably elongated, however, the 
C 2 V symmetry was still maintained. The energy of this transition structure yields 
a barrier of 0.0117 eV/atom for the (SBP)$ — ► (TP) 6 conversion reaction. The 
other 3d— >3d transformation is the (TP) 6 — ♦ (PP) 6 transition. The three tetra- 
hedra (which were the building blocks of (TP) 6 ) in this saddle point structure, 
were highly distorted. This is the only second order saddle point [9] found in this 
study. Among many paths connecting equilibrium (TP) 6 and (PP)« configurations 
no first order saddle point was found with a lower energy. The calculated energy 
value for this saddle point configuration produces a barrier of 0.0078 eV/atom for 
the (TP ) 6 — ♦ (PP) 6 transformation. The other saddle point structure, considered 
in this work, belongs to the (PP) 6 — * (CP ) 6 transition which is a 3d— >2d transfor- 
mation. This saddle point structure is a highly distorted pentagonal pyramid. The 
calculated barrier height for the (PP) 6 — ► (C'P)e transition was found to be 0.0356 
eV / atom. 

For (TBP)s — ► (SP) 5 and (SBP)$ — ► (TP)$ transformations, barrier heights 
(in both forward and reverse directions) were found to be quite low. This implies 
that isomerization reactions can easily take place at even low temperatures, and 
isomers may coexist in the gas phase for a certain length of time. Energy barriers 
for 3d— >2d transformations, for the cases of AI5 and Al6, were found to be generally 
higher than 3d — * 3d transitions, even though 2d structures are energetically more 
stable than 3d structures. This indicates that 3d— »2d transformation reactions for 
AI5 and Al 6 would proceed much slower. 

In gas-phase experiments, in general, the time necessary for a complete equi- 
libration (including all the fragmentation and structural transformation processes) 
is not a simple quantity to obtain. In a cluster generation process there are many 
factors (such as the temperature, the type of the carrier gas, dimensions of the ex- 
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pansion chamber, etc.) affecting the initial formation and structural characteristics 
of metal atom clusters. Unfortunately, present condensation theories are gener- 
ally inadequate to provide a definitive prediction about the clustering conditions of 
metal vapors in free jet experiments [11]. Based on purely statistical considerations, 
however, in the very early stages of the clustering process, the probability of a 3d 
cluster formation is higher than that of the formation of a 2d cluster. After the 
equilibration period, at relatively lower temperatures, the probability of finding a 
cluster in its lowest energy configuration is, of course, the highest. For the AI5 case, 
for example, let us compare (SP ) 5 — ► (CP) 5 and (SP) 5 — ► ( TBP) 5 transformations 
which are 3d— >2d and 3d— >3d processes, respectively. If preexponential factors of 
these processes are assumed to be approximately equal, then, based on the absolute 
rate theory (say, at T ~ 300° K) the rate for the (S P)$ — > (CP ) 5 reaction is more 
than 100 times slower than the rate of the (SP ) 5 — ► (TBP ) 5 transformation. (For 
transition energies see Table 1). This indicates that even though the low energy 
structure of AI5 is planar, the newly formed AI5 may remain 3d, fluctuating back 
and forth between (SP)s and (T BP)*, configurations a great many times before 
the complete equilibration. It has already been shown that kinetics may play a 
major role in many of the processes occurring in the gas-phase [12]. In particular, 
if the lowest energy structure of a cluster is separated by relatively high barriers 
from other low-lying configurations (like AI5 and Al« cases examined in this pa- 
per), then, the equilibration time may be long and it may become an important 
consideration in cluster beam experiments. Barrier heights between isomers may 
also play an important role in the kinetics of fragmentation which is an important 
process affecting the gas-phase composition. Among many isomeric configurations 
of a cluster perhaps one (or, a few) isomeric structure(s) may undergo to a frag- 
mentation process. Accordingly, the fragmentation depends on the concentration 
of a particular isomer in the beam. The rate of formation of that particular isomer 
(which serves as a window for the fragmentation) then, becomes an important step, 
and the energy barriers separating this isomer from other low-lying configurations, 
may control the rate of the fragmentation process in the beam. 

The model potential energy function employed in this study can satisfactorily 
reproduce highly accurate ab initio results for a variety of minimum energy con- 
figurations of small A1 clusters [10]. However, care should be exercized in a more 
quantitative interpretation of the saddle point results obtained here. This is be- 
cause, at the present time, no accurate calculations based on first principles are 
available for the barrier heights corresponding to isomeric transformation of small 
A1 clusters for comparison. 
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Table 1 


Calculated Energies and Symmetry Groups for Equilibrated 
and Transition Structures of A1 Clusters 


Cluster 

Configuration 

Energy 
(eV / atom) 

Symmetry 

Al 4 

Tetrahedral (T )4 

-1.1069 

T d 

Al 4 

[(T). - (*«» 

-1.0770 

c 2v 

Al 4 

Rhombus (R ) 4 

-1.1255 

D 2h 

Al s 

Triangular bipyramid ( TBP) 5 

-1.2469 

D 3h 

Al s 

l(TBP) s -> (5P) s ]t 


c 2v 

AI5 

Square Pyramid (SP) 5 

-1.2431 

c 4v 

Als 

[(SP) 5 - (C'/») 5 ]t 

-1.2172 

Cl 

AI5 

Close- Packed (CP) 5 

-1.2744 

C 2v 

Als 

Square Bipyramid (SBP)e 

-1.3424 

Oh 

Ale 

[( SBP) 6 - (TP) 6 ]t 

-1.3307 

C 2v 

Ai 6 

Tripyramid ( TP) 6 

-1.3480 

c 2 t> 

Ale 

{(.TP), - (PP),]< 1 

-1.3402 

C. 

"Ale 

Pentagonal Pyramid (PP)e 


C 5v 











































APPENDIX II. 


Interatomic Potential Functions with Multi-body Interactions 


During the last 5-6 years, perhaps due to the availability of supercomputers, 
simulation studies have been gaining importance at an ever increasing pace. In many 
areas today, computer simulations are becoming an integral part of many investiga- 
tive procedures and provide help in understanding various problems at atomistic 
levels. This information, in turn, has been used successfully in the interpretation of 
many experimental results. Many of these simulation techniques, which are called 
atomistic level simulations, are based on empirical model potentials describing in- 
teractions among the atoms in the system [1]. 

The concept of model potentials, generally speaking, is based on the Born- 
Oppenheimer approximation [2]. If it is assumed that in the absence of external 
forces a function $(fi , •••, r/v) exists to describe the total potential energy of a system 
of N atoms as a function of their positions, then, without any loss of generality the 
function $ can be expanded as: 

$ = 4>2 + <t>3 + f <t>n H 

where, fa, fa and <t> n denote sums of the the two-body, three- body and n-body 
interactions, respectively. In this so called many body expansion of $, it is usually 
believed that the series has a quick convergence, therefore, the higher moments 
may be neglected [2,3]. Otherwise this equation can not be employed for systems 
containing more than only a few atoms. 

In the earlier calculations, in general, the higher terms including even the three- 
body part were omitted, and the total potential energy, $, was approximated only 
by the sum of two-body interactions. This approach, which may be regarded as a 
first order approximation, not only simplified the statistical mechanical formalisms 
used in calculating various thermodynamical properties, but, more importantly, 
it enabled many earlier researchers to run simulation calculations with relatively 
smaller and less powerful computers. In most of the simulation calculations which 
are carried out considering this first order approximation, Lennard-Jones type func- 
tions were employed to mimic two-body interactions [4,5]. Despite the fact that 
those so called Lennard-Jones systems may only represent rare gases where the role 
of many-body forces is minimal, they provided a very useful understanding about 
many properties and processes in a systematic way that could not be acquired easily 
by other means [1]. Recent studies, however, have indicated that, particularly in the 
case of systems containing atoms other than those with close-shell structures, this 
first order approximation is inappropriate and produces results inconsistent with 
many experiments due to neglect of many-body interactions [6,7]. In more recent 
simulation studies, therefore, in addition to two-body interactions, three-body inter- 
actions also are being considered in the calculation of potential energies [7,8]. It has 
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been well recognized that the importance of three-body interactions increases with 
increasing covalent character of the bonding among atoms in the system. While for 
rare gases the three-body part is negligibly small, for metals and in particular for 
semiconductors their contribution turns out to be extremely significant. 

In the present work we have considered four different potential energy func- 
tions which can be expressed as separate sums of two- and three-body interactions. 
A comparative analysis was conducted for the Tersoff (T), Tersoff-Dodson (TD), 
Stillinger- Weber (SW) and Pearson- Takai-Halicioglu-Tiller (PTHT) potentials to 
determine their ability to reproduce various bulk, surface and small cluster proper- 
ties of silicon. 


Potential Functions 

1. The Tersoff potential is designed specifically for calculating energetics and struc- 
ture related properties of covalently bonded systems [9]. It can easily be written in 
sums of two- and three-body interactions. The two-body part is given by: 


i i 

and the three-body part is expressed as: 


with 


t j k 

i^j k^i,j 


r /<:(*•»* )l 

n exp{n\ 2 {rij - r ifc )) 


c + exp 

'-d-cosQji^j 


This potential uses a cut-off function given by 
fc(rij) = 1 , for Tij < R — D 

= for R-D<r i} <R + D 

fc{rij) = 0 , for Tij > R + D 

Parameters of this potential evaluated for Si systems are tabulated in Table 1. In his 
original work Tersoff [9] considered T] — 1 and Aj = 2A 2 which reduces the number 
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of parameters to be determined from 11 down to 9. Two of these, R and D, are 
used to describe the cut-off function / c (r). The parameter evaluation process was 
accomplished by fitting the potential function to reproduce the following six data: 
the cohesive energy, lattice constant and bulk modulus of bulk Si, and the cohesive 
energies of the Si dimer and the hypothetical simple cubic and face-centered-cubic 
structures as calculated by ab initio methods [10,11]. 

2. In his recent paper Dodson [12] analyzed the Tersoff potential and has con- 
cluded that the potential does a good job in reproducing many properties of silicon. 
However, he has discovered one important shortcoming. According to Dodson [12], 
the energetically most favorable structure for Si, calculated by the (T) potential 
was not a diamond cubic, but rather was a bcc structure with an energy value ap- 
proximately 1 eV/atom lower. In order to overcome this shortcoming, Dodson has 
modified the (T) potential. The new form of the potential that we will refer to as 
the Tersoff- Dodson (TD) potential, has basically the same analytical form as the 
(T) potential with an exception of the parameter tj which was taken to be equal 
to unity by Tersoff. Furthermore, Aj and A 2 have been treated as independent pa- 
rameters. Accordingly, the (TD) potential contains basically 11 parameters to be 
decided. Values of the parameters for the (TD) potential calculated by Dodson are 
also given in Table 1 along with the parameters of the (T) potential. 

3. Stillinger and Weber [13] developed a model potential comprising two- and three- 
body contributions to describe interactions in solid and liquid forms of covalently 
bonded systems. The two-body part in this potential is expressed by: 

N N 

<t>2 = s S e M r h) 

i 3 
i<j 


with, 

h( r ij) = - ( r ij)~ 9 ) * for < a 

Mr'j) = 0, for r*j > a 

where, r t * = /a, and is the distance between the i’the and j’th atoms. 

The three-body part is given by: 

N N N 

= Yi '52 e M r ij’ r ik> r jk) 

i j k 

i<j<k 
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with, 


h( r iji r iki r jk) = h ( r ij> r ik’°jik) + h (r*i,r* k ,8 ijk ) + 

where, djik represents the angle at the vertex i between the sides ry and r^. The 
function h is given by: 

h ( r ij > r ik > B jik ) = A • exp j | x ( cosOjik + § ) 2 , for r*j < a 


and 

H r ij > r ik >0jik) = 0, for r ?■ > a. 

For a specific system like Si this potential has a toted of 9 parameters to be eval- 
uated. e and a are the energy and distance scaling parameters and the rest of 
the parameters are dimensionless. The parameter a serves as a cut off function 
and determines the range of two- as well as three-body interactions. Because the 
functions cut off smoothly at r^j = a no difficulties are involved in taking proper 
derivatives for simulation studies. These parameters have been evaluated first by 
Stillinger and Weber [13] for silicon ensuring that the diamond structure is the most 
stable periodic arrangement of Si atoms at low temperature. Other criteria used in 
deciding the best values for parameters were the melting point of Si and the liquid 
Si structure. Numerical values of seven dimensionless parameters obtained by Still- 
inger and Weber [13] are: A = 7.049556277, B = 0.6022245584, p = 4.0, q = 0.0, 
a = 1.80, A = 21.0 and 7 = 1.20. The most commonly used scaling parameters are: 
e = 2.16817eV and <r = 2.095lA. 

4. Among those potential functions we included in this review, the (PTHT) poten- 
tial has the simplest analytic form [14]. Two-body interactions in this potential are 
represented by a Lennard- Jones type function: 


N N 




(ir) 


12 


J 

i<j 


ij 


2(p) 


°.)6 
ij 


( 4 ) 


where, r*j denotes the internuclear distance between the particles i and j; r Q rep- 
resents the equilibrium distance and e is the two-body energy at = r 0 . The 
three-body part in the (PTHT) potential has been approximated by an Axilrod- 
Teller type triple dipole function [15]: 


N N N 

i j k 
i<j<k 


Z( 1 + 3 CosOi Cos9j Cos9k) 

(rij • Tik-Tjk ) 3 


( 5 ) 
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where, 0i,0j,9k and rij,rik,rjk represent the angles and the sides of the triangle 
formed by the three particles i, j and k, respectively. The three-body intensity 
parameter is denoted by Z. This potential contains only 3 parameters (e, r 0 and 
Z) to be evaluated. For Si these parameters were calculated as e = 2.81 7eV, r 0 = 

o o 9 

2.2915.4 and Z — 3484 eV / A from fits to the equilibrium distance of Si 2 and 
the nearest neighbor distance of the equilibrium diamond cubic structure of Si 
[14]. Furthermore, it was ensured that the diamond cubic is the lowest energy 
configuration among the common crystalline structures and, at the same time, the 
stability condition given by {d$ / dV)T=QK = 0 is satisfied. As it is expressed 
here the (PTHT) potential does not intrinsically use a special cut-off function. 
In computations, however, summations are typically carried out beyond the 3r 0 
limit thereby including up to the sixth coordination shell around a central Si atom. 


Comparisons and Discussions 

Next, we analyzed comparatively the flexibilities of these functions and assessed 
their abilities in reproducing properties of Si for bulk, surfaces and small clusters. 

Bulk properties calculated by (T), (TD), (SW) and (PTHT) potentials for 
various structures are tabulated in Table 2. Because the equilibrium volume of 
diamond cubic Si was used in parameter evaluation procedures, all four potentials 
reproduce crystalline dimensions of bulk Si(diamond) correctly. The cohesive energy 
of diamond Si was employed in the parameterization of (T) and (TD) potentials; 
therefore, these potentials can exactly reproduce the cohesive energy of diamond Si. 
However, cohesive energy of Si(diamond) calculated by (SW) and (PTHT) poten- 
tials (with parameters given in this article) deviate about 7 and 19%, respectively, 
from the experimental value. Also given in Table 2 are cohesive energy values for 
the /3-tin structure (which is a high pressure form of Si) estimated by (T), (TD), 
(SW) and (PTHT) potentials. While (TD) potential provides the best value, the 
result obtain by (T), (SW) and (PTHT) potentials deviate about 3%, 7% and 10%, 
respectively, from the ab initio results [10,11]. Structural properties (such as the 
atomic volume and c/c ratio) are rather reasonably well predicted by both the (T) 
as well as the (PTHT) potentials. Unfortunately, values calculated by the (TD) 
potential for those structural properties were not given in the reference [12]. How- 
ever, it is anticipated that these values should be quite close to ones obtained by 
the (T) potential. The cohesive energy and the atomic volume for the (3 - tin struc- 
ture calculated by the (SW) potential were taken from Figure 1 of reference [16], 
while no value for the equilibrium c/a has been found. It appears that the (SW) 
potential is quite successful in calculating the cohesive energy of the /3-tin struc- 
ture, even though the equilibrium volume was found to be somewhat larger than 
experimental value. The (SW) potential does not, however, result in the diamond 
to /3-tin phase transformation at high pressures [17]. The cohesive energies calcu- 
lated by (T) and (TD) potentials for the hypothetical graphitic structure are in 
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fair agreement with the ab initio result obtained by Yin and Cohen [18]. Values 
obtained by the (SW) and (PTHT) potentials, on the other hand, deviate from this 
ab initio value considerably (see Table 2). While the nearest neighbor distances 
obtained by the (T) and (PTHT) potentials are in good agreement with the ab 
initio results, the (SW) potential produces a somewhat larger value for the nearest 
neighbor distance. In analyzing the cohesive energy and other structure related 
properties of the graphitic form calculated by the (T), (TD) and (SW) potentials, 
the short range nature of these functions should be kept in mind. Energy calcula- 
tions performed using these potentials for the graphitic structure, for instance, do 
not include interactions between atoms located in neighboring planes. This point 
may become an important issue if one relaxes the graphitic structure in the direc- 
tion perpendicular to its basal planes. Elastic properties of Si(diamond) at different 
temperatures have been investigated by Kluge et al., [19] using the (SW) potential. 
The maximum deviation from the experimental value found in the calculated elastic 
constants was about 40% for the C44 value. This is a very good prediction consid- 
ering that elastic properties have not been employed in the parameter evaluation 
of the (SW) function. No comparable study has been conducted for the (T) and 
(TD) potentials yet. However, both of these potentials in their parameterizations 
employ the experimental bulk modulus; therefore, elastic constants are expected to 
be reproduced within acceptable limits by these functions. On the other hand, the 
bulk modulus calculated by the (PTHT) potential is about 2.5 times larger them 
the experimental value indicating that the total potential energy curve described 
by the (PTHT) as a function of the volume is too steep around its minimum re- 
gion [20]. Melting behavior of Si using the (SW) potential has been investigated 
by various researchers [13,21-24]. In their original paper Stillinger and Weber [13] 
using a fixed volume simulation procedure calculated the melting point of Si as 
0.08 in reduced units (corresponding to T m = 2012 K). This was also reconfirmed 
later by Gawlinski and Gunton [21]. Other calculations [22,23], in general, based 
on fixed pressure procedures produced T m ss 1760 K. Takai et al. [25,26], using 
a fixed pressure Monte-Carlo method based on the (PTHT) potential, calculated 
the melting point as T m « 2000 K. Both, (SW) and (PTHT) potentials produced 
an approximately 5% volume contraction upon melting which is about half of the 
experimental value [23,25]. At the present time there has been no investigation on 
the melting behavior of Si employing (T) or (TD) potentials. 

Structure and energy values for Si(lll) and Si(100) surfaces calculated by (T), 
(TD) and (PTHT) potentials are given in Table 3, along with some corresponding 
literature values. For the Si(lll) surface, relaxation energies (with respect to the 
unrelaxed ideal surface) obtained by the (T) and (TD) functions are only slightly 
lower than the literature values. The (PTHT) potential, in this case, produced an 
energy value consistent with ab initio results. The contraction in the top interlayer 
spacing of the relaxed Si(lll) surface, was correctly predicted by all three potentials. 
The magnitude of the shrinkage is also equally well predicted by (T), (TD) and 
(PTHT) potentials. On the relaxed Si(100) surface, (T) and (PTHT) potentials 
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predict dimerization of Si atoms which is consistent with experimental observations. 
Also, dimer energies and dimer bond lengths calculated by these potentials agree 
reasonably well with the literature values. For the (TD) potential Dodson [12] in his 
paper gives no numerical values for the Si dimer on the Si(100) surface. However, 
based on Dodson comments, we expect that the dimer energy and the bond length 
predicted by the (TD) potential should be similar to the (T) potential results. 
Although the (SW) potential has been used in several computer simulation studies 
of silicon surfaces, no easily usable numerical values for the relaxation energy and 
associated structural reconstruction of the Si(lll) and Si(100) surfaces have been 
found. Only, Gawlinski and Gunton [21] in their molecular dynamics study on 
the epitaxial growth of the Si(100) surface using the (SW) potential, found clear 
evidence of the dimer formation (see also [22]). 

Optimized energies and geometries of small clusters calculated by (T), (TD), 
(SW) and (PTHT) potentials using parameters given in this report, are tabulated 
in Table 4. For the dimer case, the energy and the equilibrium distance were 
correctly reproduced by the (T) potential, because these quantities were employed 
in the evaluation of parameters. Literature values for the dimer energy and the 
equilibrium separation were reported as 3.21 eV and 2.246 A, respectively [27]. For 
(TD), (SW) and (PTHT) potentials deviations in energy values of Si 2 are 12.5%, 
-32.4% and -12.1%, while deviations in the equilibrium separation were found to 
be -2.5%, 4.6% and 2.4%, respectively. For Si 2 , accordingly, the equilibrium bond 
energy and the bond distance are best represented by the (T) potential, while the 
results obtained by the (SW) potential deviate considerably from the literature 
values. For the Si 3 case, the energetically most stable structures for the (T) and 
(SW) potentials were found to be an equilateral triangular shape while (TD) and 
(PTHT) favored linear configurations. The energy difference between these linear 
and equilateral triangular states for the (PTHT) potential was found to be 0.08 eV 
while for the others, i.e., for (T), (TD), and (SW), the energy differences are about 
1. eV or more (see Table 1). High level ab initio calculations for the Si 3 case predicts 
an isosceles triangle with an apex bond angle of 80.6° as the ground state [28]. None 
of the potentials considered in this work can provide this. However, according to 
a recent article by Koutecky and Fantucci [29], the energy difference between the 
linear and triangular forms of Si 3 would be too small to be predicted correctly. 
This is in agreement with the result obtain by the (PTHT) potential. Consistently, 
for all four cases, calculations produced interatomic distances for the triangular 
case somewhat larger than the nearest neighbor distances in the linear structures. 
While the distances calculated by (T), (TD)and (PTHT) were found to be relatively 
close to each other, (SW) results turned out to be somewhat larger. Optimized 
results for the Si 4 case are rather interesting. Each of these potential functions 
produced a different minimum energy structure. The (T) potential produced a 
three-dimensional regular tetrahedral structure as the lowest state, while the (TD) 
potential favored a linear chain structure. Planar rhombus and square structures 
were obtained as the energeticaly most stable configurations from the (PTHT) and 
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(SW) potentials, respectively (see Table 4). As in the case of Si3, the smallest 
variations in the optimized energies among different Si4 geometries was produced 
by the (PTHT) potential. There are several recent investigations [29,30] on Si4 
supporting the rhombus structure to be the most favorable one, as predicted by the 
(PTHT) potential. 

The (T) potential does a remarkably good job in calculating various energy and 
structure related properties of bulk Si in different allotropic forms. It has been also 
used to estimate several point defect energies of Si(diamond). With the exception of 
the hexagonal site interstitial energy (which was found to be almost zero), calculated 
results for defect formation have been reported to be in fair agreement with ab initio 
values [9]. In addition to surface results given in Table 3, adatom binding energies on 
the Si(lll) surface have also been calculated by the (T) potential [9]. Results seems 
to be consistent with ab initio calculations. Additionally, for the Si(lll) surface, the 
(T) potential provides the 7r-bonded chain structure of Pandey [31,32] and at the 
same time it produces the (7x7) reconstructed Si(lll) surface as the energetically 
most favorable for the Takayanagi model [9]. Even though these results produced 
by the (T) potential may be open to various criticisms, because the (T) potential 
produces a metastable Si(diamond) structure (based on Dodson’s evaluation [12]), 
the agreement obtained is quite remarkable. 

According to Dodson [12], the (TD) potential predicts correctly the diamond 
structure for Si to be the most stable structure. Undoubtedly, this is an import ant 
improvement, because the (TD) potential also produces several bulk and surface 
properties of Si with the same accuracy as the (T) potential. In our view, however, 
the (TD) potential needs further testing. In particular, values for point defects, 
adatom binding energies and binding sites have not yet been calculated. Further- 
more, as indicated above, small cluster results of these two potentials sire quite 
different. In addition to the Si2 case for which (T) and (TD) potentials produce 
somewhat different equilibrium bond energies and distances (see table 4), the (TD) 
potential values for Si3 and Si 4 cases are considerably different (both energetically 
and structurally) from the (T) potential results. Neither potential produces an S14 
configuration consistent with the ab initio result. Furthermore, the energy differ- 
ence between the linear and equilateral triangular configurations of Si3 calculated 
by both potentials are quite large [29]. Neither (T) nor (TD) potentials has yet been 
used to calculate temperature dependent quantities of Si. It will be very interesting 
to employ these potentials, for instance, for simulating the melting process. 

Cohesive energy values calculated by the (SW) potential for different forms of 
Si are generally somewhat higher than literature values (see Table 2). However, in 
addition to its ability to predict excellent melting properties, the (SW) potential 
has been shown to produce acceptable values for the pair correlation function of 
Si at elevated temperatures [13]. Also, this potential has been used successfully by 
Kluge et al., [23] to simulate the formation of amorphous silicon by rapid quench- 
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ing. According to Dodson [17] the (SW) potential produces a vacancy formation 
energy which agrees well with the ab initio result. Also, in the same paper, Dodson 
employing a Monte-Carlo procedure determined that an epitaxial growth beyond 
one-third of a monolayer coverage on the Si(lll) surface is not well represented by 
the (SW) function. Contrary to this result, however, molecular dynamics calcula- 
tions conducted by Schneider et al. [33] produced a well stacked growth and properly 
crystallized Si structures at some intermediate temperatures. In another molecular 
dynamics study based on (SW) potential, Gawlinski and Gunton [21] simulated the 
molecular beam epitaxial growth of the Si(100) surface at different temperatures. 
Also using the (SW) function Abraham and Broughton [22] analyzed the melting 
behavior of Si(lll) and Si(100) surfaces by a molecular dynamics method. Fur- 
thermore, Landmann et al. [24] successfully simulated the melt-Si(lOO) interface 
using a molecular dynamics technique based on the (SW) potential. The (SW) 
potential has also been employed by several investigators to study small clusters 
of Si. Blaisten-Barojas and Levesque [34] investigated various charged and neutral 
clusters of Si at the low temperature limit and correlated their findings with the 
observed “magic numbers”. Recently, Feuston et al., [35] using a molecular dynam- 
ics technique based on the (SW) function determined finite temperature structures 
of Si clusters and investigated the fragmentation process. The (SW) potential pro- 
duces a planar square Si 4 structure as the energetically most stable configuration 
and favors the equilateral triangular shape trimer more than 1 eV over the linear 
Sis. As indicated in Table 2, the (SW) function produced the largest deviations 
in the equilibrium energy and bond distance for the dimer among the functions 
included in this study. 

Bulk structural properties are in general well predicted by the (PTHT) poten- 
tial. Cohesive energies of different Si forms are somewhat overestimated. However, 
the (PTHT) function correctly predicts the energetically most favorable structure 
to be the diamond cubic. Despite the fact that it produces a large bulk modulus 
value, as indicated above, the (PTHT) function has been shown to provide proper 
absolute stability criteria for the diamond cubic structure [36]. Defect energies cal- 
culated by the (PTHT) potential seem to be only in a semi-quantitative agreement 
with the values reported in the literature. The vacancy and the interstitial ener- 
gies were calculated as 1.5 and 0.2 eV, respectively, which are considerably lower 
than the literature values [37]. The (PTHT) potential has been used by the same 
group in a number of surface related investigations. Computer simulation studies 
have been conducted successfully to investigate: the effect of surface stress on the 
reconstruction of the Si(lll) surface [38], Si(lll) cleavage and the (2x1) recon- 
struction process [39], long range ledge-ledge interaction on Si(lll) surfaces [40], 
adatom diffusion and adatom-ledge interactions on the Si(lll) surface [41], and 
kink site formation energies on the Si(lll) surface [42]. In these studies the level 
of success extends from quantitative to qualitative levels when compared with the 
limited number of available experimental results. In general, it is fair to say that 
structural properties of surfaces are reasonably well predicted by the (PTHT) func- 
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tion. Although the equilibrium bond energy obtained by the (PTHT) potential 
deviates about 12% from its literature value for the gas-phase Si dimer, it produces 
the bond distance of Si2 within acceptable limits. This potential, like the (TD) 
function, favors the linear Si3 configuration. However, the (PTHT) produces a very 
small energy difference between the triangular and the linear forms of Si3 which is 
consistent with the prediction by Koutecky and Fantucci based on ab initio consid- 
erations [29]. A rhombus structure for S 14 was predicted by the (PTHT) function to 
be the energetically most stable configuration that is also consistent with reported 
ab initio results [28,29]. 

As anticipated, the type of potential energy functions used in a modeling pro- 
cedure dictates many calculated properties of a system. One of the features most 
desired of a model potential function is its transferability. This means that the 
same potential function with the same parameters is able to describe various prop- 
erties of a system in different forms [9]. The evaluation of parameters of a model 
potential function is a time consuming and cumbersome job. The difficulty, in gen- 
eral, increases with the increasing number of parameters to be determined. On the 
other hand, a function with many adjustable parameters is expected to do a better 
job in reproducing different properties along with a better transferability. Perhaps, 
because of the limited number of samplings considered in the present comparison, 
potentials containing more parameters do not necessarily produce remarkably better 
results. 


Conclusions 

Structural properties of bulk Si are equally well predicted by all four poten- 
tial functions considered in this report. Cohesive energies of different bulk forms 
are best represented by (T) and (TD) functions, while calculations based on (SW) 
and (PTHT) functions exhibit some deviations. Bulk moduli of silicon estimated 
using (T), (TD) and (SW) potentials agree well with values reported in the lit- 
erature. However, the (PTHT) function produces a bulk modulus which agrees 
only qualitatively with the experiment. The melting point and the change in vol- 
ume upon melting of Si, are well predicted by both (SW) and (PTHT) functions. 
High temperature properties of Si, including the melting behavior, however, are not 
yet available for (T) and (TD) potentials. Reconstruction patterns and energetics 
of Si ( 1 1 1 ) and Si(100) surfaces were reasonably well predicted by (T), (TD) and 
(PTHT) potentials while no calculations for the same properties were found for the 
(SW) function. This conclusion on the predictive power of these functions, how- 
ever, must be treated with caution, because the present comparison is based only 
on a limited number of calculations which are available in the literature for more 
than one type of potential function. Small cluster properties (for Si3 and S14 , in 
particular) appear to be best reproduced by the (PTHT) function while the (T) 
potential provides the best result for the Si2 case. 
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Table 1 . Parameters of the Tersoff (T) and Tersoff- Dodson (TD) Potentials for Silicon [9,12]. 
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Table 2. Bulk properties of silicon calculated by (T), (TD), (SW) and (PTHT) potentials. 



(*) No calculations are available for these cases. 





































Table 3. Surface properties of Si calculated by (T), (TD) and (PTHT) potentials. 





Table 4. Small cluster properties of Si calculated by (T), (TD), (SW) and (PTHT) functions. Total binding energies 
minimum distances were denoted by £?& and r m j n , respectively. 
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